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ABSTRACT 

A  fully  Sinc-Galerkin  method  for  the  numerical  recovery  of  spatially  varying  diffusion 
coefficients  in  linear  parabolic  partial  differential  equations  is  presented.  Because  the  pa¬ 
rameter  recovery  problems  are  inherently  ill-posed,  an  output  error  criterion  in  conjunction 
with  Tikhonov  regularization  is  used  to  formulate  them  as  infinite-dimensional  minimization 
problems.  The  forward  problems  are  discretized  with  a  sine  basis  in  both  the  spatial  and 
temporal  domains  thus  yielding  an  approximate  solution  which  displays  an  exponential  con¬ 
vergence  rate  and  is  valid  on  the  infinite  time  interval.  The  minimization  problems  are  then 
solved  via  a  quasi-Newton/trust  region  algorithm.  The  L-curve  technique  for  determining 
an  appropriate  value  of  the  regularization  parameter  is  briefly  discussed,  and  numerical  ex¬ 
amples  are  given  which  demonstrate  the  applicability  of  the  method  both  for  problems  with 
noise-free  data  eis  well  eis  for  those  whose  data  contains  white  noise. 


’This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Con¬ 
tract  No.  NASl-18605  while  the  author  was  in  rooiHonrn  -t  the  Insiltui'  for  Cor.ipuiet  Applications  in 
Ci-ieiite  and  Engineering  (ICAS**)).  NASA  Langley  Research  Center,  Hampton,  VA  23665. 
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1  Introduction 


In  this  paper,  a  fully  Sinc-Galerkin  method  is  introduced  for  the  numerical  recovery  of 
variable  diffusion  coefficients  in  linear  parabolic  partial  differential  equations.  To  illustrate 
the  method,  consider  the  problem  of  estimating  the  spatially  varying  parameter  p{x)  in  the 
diffusion  equation 

u(0,<)  =  u(l,f)  =  0,  f>0  (1.1) 

u(x,  0)  =  0,  0  <  a:  <  1 

given  measurements  of  the  data  at  the  points  {(ajp,  fq)}pli’...'"p  in  (0, 1)  x  }R^ .  As  noted  in 
[1],  problems  of  this  type  arise  in  applications  ranging  from  physiological  modeling  to  sea 
sediment  analysis. 

For  many  applications,  it  is  physically  reasonable  to  assume  that  p  is  continuous  on  [0, 1] 
and  to  let  the  admissible  parameter  set  Q  be  defined  by 

Q  =  {p  €  /f‘(0, 1)  :  p{x)  >  po  >  0). 


With  this  definition,  the  existence  of  a  unique  solution  u  to  the  forward  problem  can  be 
obtained  on  any  fixed  time  interval,  (0,t],t  >  0,  for  /  sufficiently  smooth. 

The  objective  of  the  parameter  recovery  problem  is  to  choose  p‘  £  Q  so  that  the  solution 
u*  of  (1.1)  corresponding  to  p*  agrees  with  the  true  state  u.  In  general  however,  the  true 
state  u  is  not  known  and  measurements  are  taken  instead  from  an  observation  space  Z.  In 
this  paper,  the  data  are  taken  to  be  point  evaluations  and  the  observation  space  Z  is  defined 
to  be  Z  =  The  observation  operator  C  :  C'((0,  i)  X  (0,t])  — »  Z  is  then  given  by 

C0  =  {T/^(xp,f,)}’il:;;;;’.  (1.2) 


The  “idealized”  recovery  problem  may  then  be  formulated  as  follows:  determine  p  ^  Q  so 
that 


Cu(-,-,p)  =  d 
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where  d  is  used  to  denote  the  data.  Since  the  forward  problem  is  well-posed,  the  parameter 
recovery  problem  may  be  formulated  as 

KM  =  d  (1.3) 

where  the  nonlinear  operator  K  is  defined  by 

K{p)  =  ec-\p)f. 

The  problem  (1.3)  is  impractical  to  solve  for  several  reasons.  As  indicated  in  [12],  the 
problem  is  ill-posed  in  the  sense  that  solutions  p  (provided  they  exist)  may  not  depend 
continuously  on  the  data  d.  Hence,  discretized  versions  of  this  problem  are  likely  to  be  highly 
ill-conditioned.  Consequently,  some  sort  of  regularization  (i.e.,  stabilization)  is  required  to 
obtain  an  accurate  approximation  for  p. 

The  regularization  technique  that  is  used  is  Tikhonov  regularization  [19]  and  the  problem 
(1.3)  is  replaced  by  the  minimization  problem 

minT«(p)  (1.4) 

P€V 

where  the  Tikhonov  functional  is 

T.(p)s^^{\\IC{r)-df  +  aJM)- 

Here  a  >  0  is  a  regularization  parameter,  which  controls  the  tradeoff  between  goodness 
of  fit  to  the  data  and  stability.  The  penalty  functional  J{p)  provides  stability  and  allows 
the  inclusion  of  a  priori  information  about  the  true  parameter  p*.  Since  the  parameter  is 
assumed  to  be  “smooth”  in  the  sense  that  p  e  H^{0, 1),  the  penalty  functional  is  taken  to 
be  the  norm 

‘^(P)  =  llrili  =  ^  \p'{x)Yv{x)dx  -I-  \p{x)fv{x)dx.  (1.5) 

The  parameter  e  is  of  the  order  10~®  and  the  weight  v  is  taken  to  be  the  positive  function 
v{x)  =  x{\  —  x).  The  reasons  for  weighting  the  integral  as  well  as  including  the  second 
term  and  enforcing  J{p)  to  be  strictly  positive  will  be  discussed  in  the  fourth  section  of  this 
paper.  By  using  arguments  similar  to  those  in  [8]  and  [15]  and  assuming  that  K{p)  is  one  to 
one,  it  can  be  shown  that  with  this  definition  for  J{p),  the  solutions  p„  to  (1.4)  converge  as 
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the  regularization  parameter  a  — >  0  and  as  the  perturbations  in  the  data  and  operator  tend 
to  zero. 

Due  to  the  infinite  dimensionality  of  Q  and  that  of  the  state  space,  the  problem  (1.4) 
is  an  infinite-dimensional  minimization  problem.  In  order  to  develop  a  practical  numerical 
scheme,  the  problem  must  be  replaced  by  a  sequence  of  finite-dimensional  problems;  that 
is,  one  must  approximate  the  operator  K,  and  minimize  the  functional  Ta  over  a  finite¬ 
dimensional  admissible  subspace  of  Q. 

The  evaluation  of  /C(p)  requires  the  solution  of  the  partial  differential  equation  (1.1). 
Similar  PDE’s  must  be  solved  to  obtain  the  components  of  the  derivative  fC'{p).  The  con¬ 
struction  of  an  approximate  solution  to  these  forward  problems  commonly  begins  with  a 
Galerkin  discretization  of  the  spatial  variable  with  time-dependent  coefficients.  This  yields  a 
system  of  ordinary  differential  equations  which  is  solved  via  differencing  techniques.  Due  to 
stability  constraints  on  the  discrete  evolution  operator,  low-order  methods  with  small  time 
steps  are  often  required  to  obtain  accurate  approximations.  Moreover,  this  time-stepping 
must  be  repeated  at  each  step  in  the  minimization  of  (1.4).  A  final  difficulty  lies  in  the  need 
to  interpolate  at  data  points  which  do  not  coincide  with  the  nodes  of  the  ODE  solver. 

In  contrast,  the  method  of  this  work  implements  a  Galerkin  scheme  in  time  as  well  as 
space  thus  bypassing  many  of  the  difficulties  associated  with  time-stepping  methods  in  the 
context  of  inverse  problems.  This  possibility  v/as  first  explored  in  [12].  In  contrast  to  the 
methods  of  that  work  however,  both  the  spatial  and  temporal  basis  functions  are  taken  to 
be  compositions  of  sine  functions  with  suitable  conformal  maps. 

By  discretizing  the  forward  problem  in  this  manner,  the  optimal  exponential  convergence 
rate  is  exhibited  throughout  the  infinite  time  domain,  even  in  the  presence  of  boundary 
singularities.  The  validity  and  exponential  convergence  rate  of  the  approximate  solution 
throughout  all  time  is  especially  important  in  those  problems  in  which  the  data  is  sampled 
at  large  temporal  values  tg.  Furthermore,  the  sine  quadrature  rules  yield  coefficient  matrices 
which  arr  fficienlly  constructed  for  the  forward  problem  and  easily  updated  when  the  for- 
'^'»rd  techr.iq')'**'  =>’'e  Tiiuiyea  in  a  parameter  recovery  scheme.  The  efficiency  of  the  inverse 
scheme  is  further  augmented  by  the  fact  that  the  component  matrices  used  in  formulating 
the  finite-dimensional  penalty  functional  are  identical  to  those  used  when  constructing  the 
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forward  coefficient  matrices  and  hence  need  to  be  formed  only  once.  The  efficiency  and 
accuracy  of  the  forward  solver  and  the  ease  of  formulating  the  penalty  functional  are  then 
manifested  in  the  inverse  algorithm  for  a  large  class  of  problems. 

The  foundations  of  the  Sinc-Galerkin  method  and  the  fundamental  quadrature  rules  are 
described  in  Section  2.  A  thorough  review  of  sine  function  properties  can  be  found  in  [17]  and 
[18].  In  the  third  section  of  this  paper,  the  Sinc-Galerkin  system  for  the  forward  problem  is 
constructed  and  implementation  details  are  discussed.  The  section  closes  with  the  discussion 
of  a  very  robust  and  accurate  algorithm  for  solving  the  resulting  algebraic  system.  Section  4 
includes  the  finite-dimensional  minimization  problem  with  the  discussion  centering  around 
the  construction  of  the  various  components  of  the  Tikhonov  functional.  By  taking  advantage 
of  sine  function  properties,  efficient  routines  for  approximating  the  nonlinear  operator  )C(p) 
and  the  penalty  functional  S{p)  are  developed.  In  the  next  section,  a  quasi-Newton/trust 
region  scheme  is  outlined  for  solving  the  finite-dimensional  minimization  problem.  The  paper 
concludes  with  a  section  containing  numerical  examples.  Of  the  many  examples  tested,  those 
discussed  in  this  section  best  exhibit  the  features  necessary  for  the  practical  implementation 
of  the  Sinc-Galerkin  method.  A  brief  discussion  of  the  Generalized  Cross  Validation  (GCV) 
and  T-curve  techniques  for  choosing  the  regularization  parameter  a  is  given  at  the  beginning 
of  the  section,  and  the  applicability  of  these  techniques  in  conjunction  with  the  Sinc-Galerkin 
method  is  demonstrated  by  the  numerical  results.  Finally,  results  are  included  both  from 
data  sets  with  white  noise  and  from  sets  to  which  no  noise  was  added.  As  shown  in  these 
examples,  the  Sinc-Galerkin  method  works  equally  well  in  both  cases. 


2  Sine  Function  Properties 


For  the  Sinc-Galerkin  method,  the  basis  functions  are  derived  from  the  Whittaker  cardinal 
(sine)  function 


.  .  .  sm(7rx) 

sinc(x)  = - — oo  <  X  <  oo 

TTX 


and  its  translates 


S{k,  h){x)  =  sine 


h>0. 


For  h*  =  three  adjacent  members  of  this  sine  family  {S{k,  h*){x),  ^  =  — 1,0, 1)  are  shown 
in  Figure  1. 


Figure  1.  Three  Adjacent  Members  (^S{k,h*){x),k  =  — 1,0,  l,h*  =  of  the  Translated  Sine 
Family. 

To  construct  basis  functions  on  the  intervals  (0,1)  and  (0,oo),  respectively,  consider  the 
conformal  maps 

^w  =  '"(r^)  (2') 
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and 


T(iw)  =  ln{w). 


The  map  4>  carries  the  eye-shaped  region 


-  X  +  iy  : 


<d< 


onto  the  infinite  strip 

Ds  =  =  C  +  in  ••  \n\<d< 

Similarly,  the  map  T  carries  the  infinite  wedge 

TT  ^ 

Dw  =  {w  =  t  +  is  :  |ar5(ty)|  <  d  <  — } 
onto  the  strip  D5.  These  regions  are  depicted  in  Figure  2. 


Figure  2.  The  Domains  Ds,Db,  and  Dw. 
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The  sine  gridpoints  Zk  G  (0, 1)  in  De  will  be  denoted  Xk  since  they  are  real.  Similarly,  the 
gridpoints  wi,  G  (0,  oo)  in  Dy/  will  be  denoted  t*.  Both  are  inverse  images  of  the  equispaced 
grid  in  Ds]  tl  -it  is, 

gkh 

Xk  =  r\kh)  =  (2.6) 

and 


tk  =  r-\kh)  = 


(2.7) 


To  simplify  notation  throughout  the  remainder  of  this  section,  the  pairs  <f),  De  and  T,  Dw 
are  referred  to  generically  as  Xj  D.  It  is  understood  that  the  subsequent  definition,  theorems, 
and  identities  hold  in  either  setting.  Furthermore,  the  inverse  of  x  is  denoted  by  rp. 

The  important  class  of  functions  for  sine  interpolation  and  quadrature  is  denoted  B{D) 
and  defined  next. 

Definition  2.1.  Let  B{D)  be  the  class  of  functions  F  which  are  analytic  in  D,  satisfy 


f  \F(z)dz\  — >  0,  t  f  ±oo 

where  L  =  {is  :  Is]  <  d  <  and  on  the  boundary  of  D  (denoted  dD)  satisfy 

N(F)  =  I  \F(z)dz\  <  00. 

JdD 

The  following  theorem  for  functions  in  B{D)  is  found  in  [16]. 


Theorem  2.1.  Let  F  be  (0,1)  or  (0,oo)  when  x  =  ^  or  T,  respectively.  If  F  ^  B{D)  and 
Zj  =  fp{jh)  =  j  =  0,  ±1,  ±2,  •  •  • ,  then  for  h  >  0  sufficiently  small 


i 


F{z)dz-h  f; 


(2.8) 


Theorem  2.1  illustrates  the  exponential  convergence  rate  which  is  a  trademark  of  the  sine 
methods.  There  is  a  common  occasion  when  it  is  possible  to  evaluate  the  infinite  series 
appearing  in  (2.8),  namely  when  integrating  against  S{k,h)  o  x-  In  general,  however,  the 
series  must  be  truncated.  With  additional  hypotheses,  it  is  proven  in  [11]  and  [17]  that  the 
truncation  need  not  be  at  the  expense  of  the  exponential  convergence. 
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Theorem  2.2.  Assume  F  6  B{D)  and  that  there  exist  positive  constants  K,a,  and  B  such 
that 

[  T  G  ^((— oo,0)) 

<  K  I  (2.9) 

^-P\x(-r)\^  T  G  V’([0)  OO)). 

Then  for  h  sufficiently  small 


nil 

X'{r) 


L 


a 


(2.10) 


Theorems  2.1  and  2.2  are  used  to  establish  a  uniform  error  bound  when  constructing  an 
approximate  solution  to  the  forward  second-order  time-dependent  problems.  The  application 
of  these  quadrature  theorems  is  facilitated  by  the  identities 


[5(p,/i)  ox(2)] 


1 ,  i  =  p 
0,  i  p, 


(2.11) 


~S{p,h)Qx{z) 


0, 


I  (*  -  P)  ’ 


i  =  p 
i  ^  p 


(2.12) 


and 


nL 

dx^ 


S{p,h)ox{z) 


Z  =  l, 


which  denote  the  evaluation  at  the  gridpoint  Zi 
derivatives  with  respect  to  the  map  x- 


7r‘ 

~T' 


I  =  p 

,  *  ^  P 


(2.13) 


(i  -  pf 

of  the  sinc-map  compositions  and  their 
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3  The  Forward  Problem 


Consider  the  second-order  parabolic  problem 

Cu(x,i)  =  c<i<l,l>0 

u(0,t)  =  u(l,t)  =  0,  t>0  (3.1) 

u(xy  0)  =  0,  0  <  I  <  1  . 

To  define  the  Sinc-Galerkin  approximation  to  (3.1),  let  5i(a:)  =  S{i,hj;)  <p{x)  and 
S'{t)  =  S{j,ht)  0  T(<),  and  take  the  basis  to  be  ' where 

5o(x,o  =  5,(x)s;(0. 

The  approximate  solution  is  then  defined  by  way  of  the  tensor  product  expansion 

Nx  N,  +  \ 

UijSij{x,t),  (3.2) 

rnt  =  Mt  + Ni-^l  . 

The  mx  •  n!(  unknown  coefficients  {u,j}  are  determined  by  orthogonalizing  the  residual  with 
respect  to  the  set  of  sine  functions  This  yields  the  discrete  Gaierkin 

system 

{Cu^^m.-f,s,s;)  =  d  (3.31 

for  p  —  —Afr,  •  ■  ■  ,  and  g  =  —Mi,  ■  •  • ,  Nt.  The  inner  product  (•,  ■)  is  taken  to  be 

{F,G)  =  f  f  F{x,t)G{x,t)w{x,t)dxdt  (3.4) 

Jo  Jo 

with  the  weight 

w{x,t)  =  w{x)w'{t)  =  (<?i'(x))"^t(0)»  •  (3.5) 

A  thorough  discussion  motivating  this  choice  of  weight  can  be  found  in  [10]  and  [13]. 

Because  of  the  tensor  nature  of  the  approximate  solution,  the  domain  on  which  (3.1)  is 
posed,  and  the  form  of  the  inner  product,  the  discrete  system  (3.3)  can  be  formulated  by 
combining  the  discrete  systems  corresponding  to  the  one-dimensional  problems 


u{t)  =  r(/)  ,  0  <  <  <  oo 

n(0)  =  0 


(3.0) 
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{p(x)u'{x)y  =  r(x)  ,  0  <  X  <  1 

,3.7) 

u(0)  =  u(l)  =  0  . 

This  latter  approach  also  illustrates  the  sine  parameter  selections  which  are  needed  when 
implementing  the  method. 

Continuing  with  (3.6),  a  discrete  syrtem  is  formed  by  orthogonalizing  the  residual 
it„„(i)  —  r(<)  with  respect  to  Before  invoking  the  quadrature  rules,  integra- 

„ion  by  parts  is  used  to  transfer  the  differentiation  of  u  onto  S*Vt,  where  again,  iv*  = 
denotes  the  temporal  inner  product  weight.  To  guarantee  that  the  boundary  terms  vanish. 


it  is  assumed  tha* 


..  ^(0  1-  w(0 

lim  —7^  =  hm  — 7=-  =  0. 

y/i  t-*00  yjl 


Finally,  the  resulting  integrals  are  evaluated  via  (2.10)  o'  (2.8)  when  possible.  With  respert 
to  (2.9),  the  condition 


'^(0t/t(0l  < 


t\  <€(0,1) 
<“*,  <  €  [1, 00) 


guarantees  the  lioundcdness  necessary  to  truncate  the  infinite  quadrature  rule.  Witli  7  and 
6  specified  and  Mt  chosen,  the  pa  arneter  selections 


A',  =  +  (•^■8) 

where  [•]  denotes  the  greatest  integer  function,  balance  the  asymptotic  quadrature  errors  in 
(2.10)  to  at  least  ^ 

In  many  parabolic  systeins,  it  is  reasonable  to  assume  that  the  solution  decays  exoonen- 
tially  at  infinity,  that  is  that  the  solution  satisfies 


i\  <f(0,i) 
/f[l,oo) 
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or,  more  succinctly, 


\u{t)\  < 


(3.9) 


Under  this  supposition,  Lund  [11]  shows  that  the  condition  (3.8)  can  be  replaced  by 

‘  (3.10) 

The  selection  ;V(  in  (3.10)  significantly  reduces  the  size  of  the  discrete  system  with  no  loss 
of  accuracy.  It  is  also  noted  that  the  size  of  the  discrete  system  and  the  expected  error  are 
dictated  by  the  asymptotic  behavior  of  u  at  the  endpoints. 

The  discrete  system  for  (3.6)  can  then  be  formulated  as  follows.  Let  /<'),  £  =  0, 1  denote 
the  7ni  X  nit  matrices  whose  qj-th  entry  is  from  (2.11)  and  (2.12)  and  let  T>(t/)  be  the  di¬ 
agonal  matrix  with  entries  ■  •  •  ,T)(ti\/,).  The  vector  of  unknowns  u  = 

is  then  related  to  the  known  vector  f  =  [r(<_Ar, ),•••,  r(<A^,)]^  by 

Atu  =  Vi{r)-^)f  (3.11) 

where 

A  =  +  (3.12) 

Further  details  concerning  the  derivation  of  the  system  (3.11)  can  be  found  in  [10]  and  a 
tho'ough  analysis  of  the  spectrum  of  At  is  given  in  [3]. 

The  preceding  discussion  applied  to  the  problem  (3.7)  follows  a  similar  development.  The 
map  T  of  (2.2)  is  replaced  by  the  map  4>  of  (2.1)  (since  (3.7)  is  posed  in  the  interval  (0,1)) 
and  ht  is  replaced  by  h^.  Orthogonalizing  the  residual  and  two  integrations  by  parts  yields 
the  system 


/  p(j)f5p(j)  ^ - I  dx  =  [  r{x)Sp(x)  (3.13) 

for  p  =  —Mr,  •  ■  •  ,  A^i  do  guarantee  that  the  boundary  terms  vanish,  it  is  assumed  that 


(x)[  =  (^) 


(i)  =0. 


In  anticipation  of  the  parameter  recovery  problem  which  motivates  this  analysis,  the  term 
p(x)  in  (3.13)  is  expanded  as  a  linear  combination  of  sine  functions  with  two  Ilermite  like 


II 


algebraic  terms  added  to  accommodate  the  nonzero  values  of  p  at  x  =  0  and  x  =  1.  The 
finite-dimensional  approximation  of  p  then  takes  the  form 


PmA^)  =  -  x)  +  Cat^X  +  <^kSk{x) 

k=-Mi  +  l 
Nx 

=  2  . 
k--Mx: 


(3,14) 


In  the  forv/ard  problem,  the  coefficients  are  known  whereas  in  the  corresponding 

parameter  recovery  problem,  they  are  unknown  and  are  determined  via  methods  to  be  dis¬ 
cussed  in  Section  4.  The  number  of  basis  functions  used  in  the  expansion  is  chosen  so  as  to 
guarantee  a  square  coefficient  matrix.  This  is  done  to  simplify  the  implementation  of  the 
method  when  applied  to  the  PDE  (3.1)  of  interest. 

The  expansion  (3.14)  is  substituted  into  (3.13)  and  the  resulting  integrals  are  evaluated 
via  (2.10)  or  (2.8)  when  possible.  As  shown  in  [13],  the  decay  condition  (2.9)  equates  to  the 


condition 


lu(x)P(x)|  <  L 


^  €  (0,  \) 


(1  -  x)^+2,  X  G  [i  1) 


where 


F(x)  =  p(x)  -  p(0)(l  -  x)  -  p(l)x. 


This  may  be  replaced  by  the  more  stringent  condition 

|u(x)P(x)|  <  /fi"'*'j(l  — 


(3.15) 


The  asymptotic  errors  are  then  balanced  by  the  choices 


N.  =  I^M.  +  1| 

where  [•]  again  denotes  the  greatest  integer  function.  Note  that  if  is  an  integer,  this 
integer  can  be  selected  for  Nf 
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With  u,f,  and  (  =  0, 1,2,  defined  as  before,  the  system  for  (3.6)  may  be  written  as 


A{p)it  = 'D((<f>')  i)r 


(3.16) 


where 


Mp)  --  -i/w  + 


(3.17) 


The  notation  T>[p^)  and  VOp^i)  denotes  the  diagonal  matrices  containing  the  components  of 
the  vectors  and  which  are  defined  as  follows.  First 

p^.  =  'i/c 

where  c  =  [c_Afj,  •  •  • ,  c/v,]^  and  ^  has  the  block  structure 

^  =  :  0R]m,X7n, 


=  [(1  -  i-A/,),  •••,(!  -  ar.Vx)]^ 

and 

i’R  =  [x_Af,,-  •  • 

Again,  the  tUx  x  (m^  —  2)  matrix  has  components  as  defined  in  (2.11)  with 

—  Mx  <  <7  <  Aj  and  —Mx  +  1  <  J  <  iVj.  —  1.  Also, 

p^>  =  '^'c 


where 


Here  T  =  [1,  •  •  • ,  1]^,  is  mx  x  nix,  and  is  nix  x  {nix  —  2)  with  components  as 

defined  in  (2.12). 

As  shown  in  [13],  the  system  (3.16)  yields  an  approximate  solution  which  exhibits  expo¬ 
nential  convergence  to  the  solution  u  of  (3.7).  Further  details  concerning  the  derivation  of 
the  .system  as  well  as  additional  quadrature  hypotheses  can  be  found  in  this  reference. 
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The  above  results  for  the  one-dimensional  problems  (3.6)  and  (3.7)  can  then  be  pieced 
together  to  form  the  Sinc-Galerkin  system  for  the  time-dependent  parabolic  problem  (3.1). 
The  resulting  discrete  system  is  built  from  the  matrices  At  (in  (3.12))  and  A(p)  (in  (3.17)) 
of  the  one-dimensional  problems.  The  parameter  selections  are  still  necessary  and  all  that 
remains  is  to  asymptotically  balance  the  resulting  errors  from  each  one-dimensional  problem. 

When  the  decay  conditions  (3.9)  and  (3.15)  are  combined  to  yield 

\P{x)u{x,t)\  <  /Cx“+^(1  -  (3.18) 


then  the  choices 


nd 


—  O  ,4  > 

aMx 
ht  —  hxy 


Mt  = 


a 


7 


Mx  +  1 


) 


and 

for  the  stepsizes  and  summation  limits  balance  the  asymptotic  errors.  If  one  takes  d  ~  then 
the  above  choices  yield  an  asymptotic  error  rate  of  order  for  the  quadratures. 

Given  Mx,  Nx,  Mt,  Nt,  and  h  =  hx  =  ht  &s  defined  above,  the  discrete  system  for  (3.1)  is 


A{p)VV  ((t)-^)  -1-  V  ((<^')'^)  =  G  (3.19) 


where 

G  =  D((f)-S)fP  ((+)-*). 

The  diagonal  matrices  V  ^(C^) ”*)  have  sizes  x  rrix  and  mt  x  m,,  respec¬ 

tively.  The  nij.  X  rUf  matrices  U  and  F  contain  the  unknowns  {ujj}  and  the  known  values 

The  discrete  Sinc-Galerkin  system  (3.19)  can  then  be  .solved  for  U  via  a  generalized  Schur 
decomposition  (page  396  of  [6]).  As  guaranteed  by  the  results  of  Moler  and  Stewart  [14], 
there  exist  unitary  matrices  Q\,Zi,Q2,  and  Z2  such  that 
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Q\A{p)Z,  =  P 
Q\V  =  R 

g;p((t)-5)z2  =  5 

Q;AtZ,  =  T 


where  P,R,S,  and  T  are  upper  triangular.  If  F  =  Z^UZi  and  C  =  Q\GQ2-,  then  (3.19) 
transforms  to 


PFT*  +  RYS'  =  C. 


By  comparing  the  ^-th  columns,  one  finds  that 


j=k  j=k 


which  yields 

{hkP  +  s,kR)yk  =  ^k  -  P  E  hjyj-R  E  s^^yj  (3.20) 

j=fc+i  i=fc+i 

(for  convenience,  it  is  assumed  that  all  matrices  are  n  x  n  and  indexed  from  1  to  n).  With 
the  assumption  that  the  matrix  {tkkP  ^kkR)  is  nonsingular,  the  solution  to  (3.20)  is  easily 
found  by  recursively  solving  triangular  systems. 

Although  this  algorithm  does  require  complex  algebra,  it  is  quite  efficient  and  requires  no 
assumptions  concerning  the  diagonalizability  of  the  component  matrices.  It  should  be  noted 
that  a  “real”  version  of  this  algorithm  also  exists  [5].  In  this  latter  algorithm,  Qi,  Zi,Q2,  and 
Z2  are  orthogonal  with  P,S  queisi-upper  triangular  and  R,T  upper  triangular.  As  proven 
in  [5],  the  real  algorithm  is  extremely  stable  and  numerical  tests  have  indicated  that  the 
complex  algorithm  described  above  is  also  robust. 
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4  The  Finite-Dimensional  Minimization  Problem 


As  noted  in  the  introduction,  the  minimization  problem 


mirT.(p) 


where 


'To{p)  =  ^{\\IC{p)~df  +  a\\prQ}  (4.1) 

is  infinite-dimensional  and  hence  must  be  replaced  by  a  sequence  of  finite-dimensional  prob¬ 
lems  before  a  viable  numerical  scheme  may  be  developed.  Following  from  (3.14),  the  ap¬ 
proximating  admissible  parameter  sets  are  taken  to  be 

Qmi  =  I  Pmi  •  Pmi{x)  =  <^A:Ca:(^)  / 

I  k=~M.  ) 

where 

( 

1  —  X,  k  =  —Mx 

<  5'fc(x),  -Mx  -\-l<k<Nx~l  (4-2) 

X,  k  —  Nx 

and  5fc(x)  =  S{k,hx)  o  (f>{x).  'I'he  associated  finite-dimensional  optimization  problem  can 
then  be  formulated  as 


<k{^)  = 


min  fa{Pm,) 

Pmx  €  Vmx 


(4.3) 


for 


\  -  d  11=  +  a||p„J|J,}  .  (4.4) 

The  approximation  K{pmJ  ■  to  fC{p)  is  obtained  by  applying  the  point  eval¬ 

uation  operator  C  in  (1.2)  to  Um^m,  in  (3.2).  If  the  set  of  observation  points  {(a'p,  <<,)}p=l’...’^p 
can  be  represented  as  a  tensor  product  of  spatial  and  temporal  points,  then  /v(pm,)  has  the 
repre.sentation 

K(p^.)  =  Cc^U)  (4..'i) 

where  the  matrix  U  solves  (3.19).  'ITie  matrix  concatenation  co{U)  is  the  vector  in 
which  is  obtained  by  successively  stacking  the  columns  of  the  Wx  x  matrix  U .  C  is  an 
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(rip-nq)  X  {rrix-mt)  evaluation  matrix  which  can  be  formulated  follows.  Define  the  rip  x  rux 
spatial  evaluation  matrix  Ex  to  have  components 


[Ex]p,i  =  Si{xp),  1  <  p  <  np,  -Mx  <  i  <  Nx 


and  define  the  n,  x  tt?  ,  temporal  evaluation  matrix  Et  to  have  components 

[£;<],, j  =  l<?<n„  -Mt<j<Nt. 

Then 

C  =  Et®  Ex. 

It  is  noted  that  if  the  set  of  observation  points  is  not  rectangular  as  described  above,  then 
point  evaluation  can  be  done  directly  via  (3.2).  This  latter  option  is  less  efficient  however, 
than  that  defined  in  (4.5). 

The  discrete  penalty  functional  |lPm*||g  is  formed  by  substituting  the  expansion  (3.14) 
into  the  definition  (1.5).  This  yields 

llPmJIg  [p'^^{x)]Mx)dx  +  e  [p„,(x)]^t;(x)<ix  w  c^Qc 

where  the  nix  x  matrix  Q  —  Qd  +  Qj  has  components 

[Qd]kt  «  /  Cfc(2;)C/(a^)^(a;)dx,  -Mx  <  k,£  <  Nx 

Jo 

and 

[Qf]ke  «  £  /  C*(a^)C<(a:)i’(x)dx,  -Mx  <  k,£  <  Nx- 

Jo 

The  matrix  entries  are  approximations  in  the  sense  that  sine  quadrature  is  used  to  evaluate 
many  of  the  integrals. 

For  the  choice  of  basis  functions  in  (4.2),  the  matrix  Qd  is  given  by 

Qd  = 


<}d 


<id 

Qd 

-qd'^ 


-Rd 

1 

6 
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Integration  by  parts  and  the  application  of  the  sine  quadrature  formula  (2.10)  yields  the 
{rris  —  2)  X  (mj.  —  2)  matrix 


0.  = 

where  again  denotes  the  matrix  whose  ^j-th  entry  is  from  (2.13).  The  zer  ung  of  all 

other  quadrature  terms  is  a  result  of  the  choice  i>(i)  =  =  x(l  —  x).  The  {rux  —  2)  x  1 

vector  qd  has  components 


=  hx{xk  -  3Xfc  +  2x^),  -Mx  +  l<A:<iVx-l 


and  is  again  obtained  via  sine  quadrature.  Because  is  negative  definite  (see  [16]),  the 
matrix  Qd  is  nonnegative  definite.  The  zero  eigenvalue  results  from  the  fact  that  the  first 
and  last  columns  of  Qd  differ  only  in  sign. 

Direct  integration  and  sine  quadrature  are  also  used  to  obtain  the  matrix 


Here 


Qf  =  e 


1  -♦  T  1 

20  9/1  30 

9/1  Qs  9fT 

.JL  s,  T  ± 

30  v/r  20 


Qj  =  hxV{x\l-xy) 


where  D(r/)  again  denotes  the  diagonal  matrix  with  entries  77(x_Mx)>  •  •  • ,  The  vectors 

qji  and  qjr  have  components 

[qfi]k  =  hxxl{l  -  Xkf 

and 

[9/r]lc  =  hxXk{l  —  Xfc) 

for  k  =  ~Mx  +  1,  •  •  • ,  iVx  —  1.  The  matrix  Q/  is  strictly  positive  definite. 

Although  the  matrix  Q  is  full,  it  is  very  efficient  to  construct  since  the  Toeplitz  matrix 
is  also  needed  in  the  forward  solver.  For  e  >  0,  Q  is  symmetric  and  positive  definite  and 
hence  has  a  Cholesky  decomposition  Q  =  R  where  R  is  upper  triangular.  It  then  follows 
that  the  penalty  term  llPm^llg  yields  the  quadratic  form 

c'^R^Rc=^  |(/2c  II*  (4.6) 


I 
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where  j|  •  ||  denotes  the  Euclidean  norm.  As  will  be  shown  in  the  next  section,  this  factor¬ 
ization  admits  a  particularly  useful  diagonalization  of  the  corresponding  finite-dimensional 
minimization  problem.  It  also  facilitates  the  plotting  of  the  L-curve  to  determine  a  suitable 
regularization  parameter  a  (see  Section  6  and  [7]). 


5  The  Trust  Region  Scheme 

In  the  discussion  of  this  section,  it  is  useful  to  highlight  the  dependence  of  the  operators  in 
(4.4)  on  the  unknown  vector  c  =  [c_M,,-  •  •  (see  (3.14)).  Letting 

K{c)=^k{prr.Ac))-=^Cc^U) 


and  noting  (4.6),  the  optimization  problem  (4.3)  may  be  replaced  by 


min  Taic) 


(5.1) 


where 


1 


T^{c)^^{\\Kic)-dr  +  a\\Rcr}. 


To  obtain  a  minimizer  for  the  nonlinear  functional  Ta,  a  queisi-Newton/trust  region  scheme 
is  used  (see  [2]). 

The  basis  for  this  approach  is  the  iteration 

Ck+l  =  Ck  +  Sk 


where  s*  solves  the  constrained  minimization  problem 


i{||  A-(c-i)  +  I<\c,)h  -  d  f  +  a||7?(c-»  +  j’Of) 


(5.2) 
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subject  to  |l/?5*fc||  <  <5jt.  The  trust  region  radius  6k  is  chosen  so  that  the  quadratic  model 
adequately  reflects  the  behavior  of  Tq  within  the  trust  region;  that  is,  6k  is  chosen  so  that 

la{ck  +  Sk)  «  ^-{\\K{^k)  +  I<'{^k)sk  -  d  11^  +  a\\R{ck  +  s,)in 

whenever  <  6k.  The  minimization  problem  (5.2)  is  solved  using  an  approach  similar 

to  that  in  [9].  The  problem  is  first  diagonalized  using  the  Singular  Value  Decomposition 
(SVD).  With  the  change  of  variables 

3  =  Rsk, 

the  objective  functional  in  (5.2)  becomes 

i{||/t5-tiP+a||c  +  sf} 

where  A  =  K'[ck)R~^ ,  h  —  d  —  K{ck)  and  c  =  Rck.  Let  A  have  the  SVD 

/I  =  UDV^ 


where  n,),  are  orthogonal  and 


[Di 


npn^)xm,Ji,j 


0, 


if  i  =  j  and  i  <  rux 
otherwise. 


Here  denotes  a  singular  value  of  A.  The  second  change  of  variables 


s  =  V^s,  b  =  U^b,  c=V 


yields  the  diagonalized  problem 

+  "11^  +  (5-^) 

subject  to  ||s||  <  6k. 

The  theory  of  constrained  optimization  is  used  to  solve  (5.3).  By  the  Kuhn-Tucker 
criterion  [4],  there  exists  a  Lagrange  multiplier  /x  >  0  such  that 

[Ds  —  6)  +  a(c  +  s)  +  /xs  =  0.  (5.4) 
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From  (5.4)  it  follows  that  (5.3)  heis  a  unique  solution  of  the  form 


5  =  s(/i)  =  {D^ D  +  (a  +  fi)l)  ^{D^b  —  ac). 

If  ||i(0)||  <  6k,  then  the  constraint  in  (5.3)  is  not  active  and  s  =  jR~*Ks(0)  solves  (5.2); 
otherwise  the  constraint  is  active  and  the  solution  to  (5.2)  is  given  by  s  =  R~^Vs{n)  where 
>  0  is  the  unique  solution  to 

M  =  P(/^)ll  -  ^Jt  =  0.  (5.5) 


An  approximate  solution  to  the  scalar  equation  (5.5)  is  then  determined  via  the  hook  step 
algorithm  in  [2]  (see  page  134).  This  algorithm  requires  both  and  As  shown  in 

[9],  the  function  g{fi)  can  be  expanded  as 


9i^i)  = 


( -  QCj  y 

V<T?  +  a  +  /i; 


—  6k 


(5.6) 


t  A 

when  bj  and  cj  are  components  of  6  and  c,  respectively.  The  derivatives  g'{fi)  are  easily 
obtained  from  the  form  (5.6). 

The  trust  region  radius  6k  in  (5.3)  is  chosen  so  that  ^^(c  )  has  sufficient  decrease  at  each 
iteration  so  as  to  guarantee  convergence  to  a  local  minimizer  of  Tq.  This  is  accomplished  via 
the  updating  algorithm  in  [2]  (page  143)  with  the  decay  requirement  taken  to  be 


Ta{ck  +  Sk)  <  r^Ck)  +  eVTaickfsk 


with  e  =  10“"*. 

An  important  numerical  issue  in  the  implementation  of  the  trust  region  scheme  is  the 
formulation  of  the  derivation  of  the  operator  K.  Here  the  derivative,  or  Jacobian,  is  an 
(np  •  n,)  X  nix  matrix  whose  i/-th  column  is  given  by 

I)f'(c  )1.  =  Jim  ilif(c+  Te.)  -  K(S)1 


where  the  standard  unit  vector  e^,  has  components 


=  6„k  = 


if  k  =  i/,  —Mr  <  k  <  Nx 
otherwise. 
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In  the  examples  that  are  presented  in  Section  G,  the  Jacobians  were  calculated  with  a 
standard  forward  difference  scheme.  This  scheme  is  easy  to  implement  and  accurate  enough 
for  the  purposes  of  the  method.  If  further  efliciency  is  desired,  a  directional  derivative 
scheme  such  as  that  described  in  [12]  can  be  used.  For  this  method,  the  trade-off  for  the 
added  efficiency  is  an  algorithm  which  is  more  difficult  to  implement  and  a  slight  loss  of 
accuracy  in  some  cases. 


6  Implementation  and  Numerical  Examples 

The  four  examples  reported  in  this  section  were  selected  from  a  large  collection  of  problems 
to  which  the  Sinc-Galerkin  method  was  applied.  The  results  are  representative  of  those 
obtained  on  other  sample  problems. 

The  first  example  demonstrates  the  application  of  the  Sinc-Galerkin  method  to  a  model 
problem  in  which  the  state  solution  was  sampled  directly;  that  is,  no  external  noise  was  added 
to  the  data.  To  demonstrate  the  feasibility  of  the  method  for  problems  with  noisy  data, 
the  same  problem  is  revisited  in  Example  5.2  but  with  pseudo-random  white  noise  added 
to  the  data.  In  Example  5.3,  the  parameter  to  be  recovered  has  a  logarithmic  boundary 
singularity  at  a;  =  0  while  the  parameter  in  Example  5.4  is  the  shifted  Gaussian  function 
that  was  considered  in  [12].  In  all  four  examples,  the  dynamics  of  the  problem  are  assumed 
to  be  modeled  by  (1.1)  with  the  forcing  function  /(x,f)  consistent  with  the  state  solution 
u{x,i)  =  x(l  -x)sin(47rx)f’c“*  and  the  true  diffusion  parameter  p.  In  each  case,  cf  =  |  (sec 
(2,3)  and  (2,5)). 
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The  errors  for  the  method  are  reported  on  tlic  set  of  uniform  gridpoints 


U  —  {^'o,  •  •  • ,  2^100} 

Zk  =  =  U,  1,- '  • ,  100 


(C.l) 


with  stepsize  £  =  With  p  and  Pm,  denoting  the  true  and  approximate  parameters 
respectively,  the  errors  are  reported  as 


The  error  and  convergence  results  arc  tabulated  in  the  form  .aaa  —  7  which  represents 
.aaa  x  10“'’'.  All  problems  were  run  with  sixteen  place  accuracy  on  a  Vax  8550. 

A  very  important  practical  consideration  is  the  choice  of  the  regularization  parameter  a 
for  a  given  (error  contaminated)  data  set.  One  would  like  to  choose  a  so  that  \\p  —  pa||  is 
minimized,  where  pa  denotes  the  a-dependent  unknown  diffusion  coefficient.  If  the  error  in 
the  data  is  random,  then  under  certain  conditions  (see  [20])  the  method  of  Generalized  Cross 
Validation  (GCV)  yields  a  statistical  estimate  of  the  size  of  |l/C(p)  —  AC(pa)l!  which  is  related 
to  Up  -  Pall .  For  Tikhonov  regularization,  this  estimate  is  given  by  the  GCV  functional 

l||ji-(g„)  -  S\V 

n-rrij,  ^  a 

.  "  + 

where  Ca  solves  (5.1).  Here  n  =  rig  -Up  denotes  the  number  of  data  points  and  {cr,}  are  the 
singular  values  of  the  operator  K\ca)R~^.  To  approximate  the  value  of  a  which  minimizes 
Up  —  Pall,  one  attempts  to  solve  the  minimization  problem 


minV'(a). 

a>0  '  ' 

Because  the  GCV  method  requires  the  singular  values  of  /{"'(ca),  it  is  relatively  expensive 
to  implement  when  m*  and  n  are  large.  A  second  disadvantage  to  this  method  for  choosing 
the  regularization  parameter  is  that  the  GCV  plots  are  often  very  flat  making  it  difficult  to 
determine  a  minimum  value  of  V{a)  and  hence  an  optimal  value  of  a  (see  Figure  6  in  the 
next  section).  Finally,  one  often  has  optimization  settings  in  which  the  GCV  hypotheses  are 
not  satisfied. 
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A  second  method  for  determining  the  regularization  parameter  is  to  plot  the  norm  of 
the  penalty  functional,  ||/fca||,  versus  the  norm  of  the  residual,  ||A' (cq,)  —  d  ||  (see  ["]).  In 
this  way  one  can  qualitatively  get  an  idea  of  the  compromise  between  the  minimization  of 
these  two  quantities.  The  scheme  for  determining  the  “optimal”  regularization  parameter 
consists  of  finding  those  values  of  a  such  that  (||/t'(cQ)  —  d  ||,  ||i?CQ|j)  lies  in  the  “corner” 
of  the  resulting  curve,  known  as  the  L-curve.  This  method  for  choosing  the  regularization 
parameter  a  is  easy  to  apply  and  often  gives  more  conclusive  results  than  ihe  GCV  method. 
Both  methods  are  illustrated  in  the  examples. 

In  all  four  examples,  the  x  1  initial  vector  ^  =  [.5,  .5,  ■  •  • ,  .5,  .5]^  was  used.  Several 
other  positive  startup  vectors  were  also  tried  with  similar  results  in  each  case.  Hence  the 
metliod  seems  to  be  quite  robust  with  respect  to  the  choice  of  the  initial  vector. 

Finally,  in  the  examples  the  symbol  a  is  used  to  denote  both  the  regularization  parameter 
(see  (5.1))  and  the  sine  decay  parameter  (see  (3.18)).  The  use  of  this  symbol  for  both 
quantities  is  well  established  in  the  literature  and  thus  difficult  to  avoid  in  this  setting.  It 
should  be  obvious  from  the  context,  however,  which  quantity  is  being  discussed  and  there 
should  be  no  ambiguity  resulting  from  the  dual  use  of  this  symbol. 


Example  6.1  In  this  example,  the  true  diffusion  parameter  is  taken  to  be  the  analytic 
function  p[t)  =  1  +  sin(7rx).  Since  the  state  solution  is  u{x,t)  —  a;(l  —  x)  sin(47r2:)f^e“',  the 
decay  condition  (3.18)  yields  the  choices  a  =  ^  =  2,  7  =  |,  and  <5=1.  Tlie  data  was  sampled 
on  a  regular  grid  C  (0,1)  x  (0,2].  Nine  equally  spaced  points  Xp  =  pAx,Ax  =  .1, 

were  taken  in  space  and  four  equally  spaced  temporal  points  tq  =  qAt,  At  =  .5,  were 
taken  for  a  total  of  n  =  54  data  points.  No  noise  was  added  so  the  data  consisted  of  direct 
measurements  of  the  state  solution.  For  varying  values  of  the  regularization  parameter  a,  the 
L-curve  is  plotted  in  Figure  3.  Note  that  the  values  a  =  10“^  through  a  =  10~''  yield  points 
(||A'(ra)  —  d  ||,||Acall)  “corner”  of  the  curve.  The  uniform  errors  for  a  =  10“®  are 

reported  in  Table  1  with  the  first  four  columns  indicating  the  index  limits  for  the  expansion 
of  the  stale  variable  and  fifth  column  indicating  the  number  of  basis  functions  used  in  the 
expansion  of  ■  The  convergence  of  the  method  is  demonstrated  both  by  the  results  in 
the  last  column  of  Table  1  and  by  Figure  4  which  shows  the  true  and  approximate  diffusion 
parameters  with,  a  =  10~®. 


Mx 

Nx 

Mt 

Nt 

rux 

I|P.(0II 

8 

8 

10 

4 

17 

.7414-0 

16 

16 

21 

7 

33 

7648  -  1 

24 

24 

39 

9 

49 

.3111  -  1 

Table  1.  Errors  on  the  Uniform  Grid  U  with  q  =  10  ®  in  Example  6.1. 
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Example  6.2  Here  the  true  parameter  and  state  solution  are  the  same  as  those  in 
Example  6.1,  and  hence  p(x)  =  14-  sin(7rx)  and  u(x,t)  =  x(l  —  x)  sin(47rx)<^e~‘.  The  same 
observation  points  were  used  but  to  this  data  however,  we  added  a  pseudo  random  noise 
vector  e  from  a  Gaussian  distribution  with  mean  0  and  standard  deviation  cr  chosen  so  that 
the  noise-to-signal  ratio  o’/||d||  =  0.001  (noise  -  0.1%  of  the  signal).  The  T-curve  and  GCV 
curves  are  plotted  in  Figures  5  and  6,  respectively.  Note  that  the  values  q  =  10“®  through 
a  =  5  X  10“®  yield  points  (||/f(cQ)  —  d*  ||,[[/Zco||)  in  the  “corner”  of  the  T-curve  whereas 
all  values  of  a  less  than  10~®  yield  apparent  minima  of  the  GCV  curve.  For  —  16,  the 
uniform  errors  obtained  with  a  =  10“^,  a  =  10“®,  and  ct  =  10~®  are  reported  in  Table  2. 
Corresponding  plots  of  the  true  and  approximate  parameters  are  shown  in  Figure  7.  Note 
that  the  “corner”  value  a  =  10“®  provides  a  good  choice  for  the  regularization  parameter 
whereas  a  =  10“®  is  not  large  enough  to  damp  out  the  contribution  due  to  the  smaller  sin¬ 
gular  values.  This  latter  observation  can  be  predicted  from  the  L-curve  but  less  easily  from 
the  GCV  plot.  Finally,  the  choice  a  =  10"^  causes  too  much  smoothing  and  information 
about  the  parameter  is  lost.  By  comparing  the  results  in  Tables  1  and  2,  it  can  be  seen  that 
the  error  in  this  example  with  a  =  10"®  and  Mx  =  16  is  virtually  the  same  as  the  error  in 
Example  6.1  with  a  =  10”®  and  Mx  =  16.  The  results  from  this  example  (iemonstrate  the 
viability  of  the  method  for  problems  with  noisy  data. 


\ 

a  =  10"®  Q  =  10"®  Q  =  10"® 

\\pum 

.2658  -  0  .7737  -  1  .4357  -  0 

Table  2.  Errors  on  the  Uniform  Grid  U  with  Mx  =  16  in  Example  6.2. 
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2.2 


Figure  7.  True  and  Approximate  Diffusion  Parameters  for  Example  6.2  with  Mx  =  16 
- (a  =  10-3), - (a  =  10-3),  •  •  •  (a  =  lO"®), - (True). 


Example  6.3  The  true  parameter  in  this  example  is  p(x)  =  1  +  |a:  +  xln{x)  which  has  a  log¬ 
arithmic  singularity  at  x  =  0.  As  before,  the  state  solution  is  u(x,t)  =  x(l  —  x)  sin(47rx){^e-‘ 
and  thus  the  decay  parameters  a  =  /?  =  2,7  =  |,  and  ^  =  1  are  consistent  with  the  condition 
(3.18).  To  demonstrate  the  method  for  another  set  of  observation  points,  nineteen  equally 
spaced  points  Xp  =  pAx,  Ax  =  .05,  were  taken  in  space  and  four  equally  spaced  temporal 
points  tq  —  qAt,  At  =  .5  were  taken  for  a  total  of  n  =  72  data  points.  No  noise  was  added  so 
the  data  consisted  of  direct  mezisurements  of  the  state  solution.  Since  the  L-curve  wais  nearly 
identical  to  that  of  Example  6.1,  the  regularization  parameter  was  taken  to  be  q  =  10“®. 
The  uniform  errors  for  this  choice  are  reported  in  Table  3  and  the  true  and  approximate  pa¬ 
rameters  are  shown  in  Figure  8.  Both  the  table  and  the  figure  demonstrate  the  convergence 
of  the  method  in  spite  of  the  logarithmic  singularity  in  the  diffusion  parameter. 
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N. 

Mt 

Nt 

m* 

1Ip4^)II 

8 

8 

10 

4 

17 

1.2171-0 

16 

16 

21 

7 

33 

0.1330  -0 

24 

24 

39 

9 

49 

0.7285  -  1 

Table  3.  Errors  on  the  Uniform  Grid  U  with  o;  =  10  ®  in  Example  6.3. 


Figure  8.  True  and  Approximate  Diffusion  Parameters  for  Example  6.3  with  a  = 
■  -  •  (M*  =  16), - (M,  =  24), - (True). 
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Example  6.4  In  this  example,  the  parameter  to  be  recovered  is  the  shifted  Gaussian  func¬ 
tion  p{x)  =  1  -f  When  combined  with  the  state  solution,  this  dictates  the 

choices  a  =  /?  =  7  =  |,  and  6  =  1  for  the  sine  decay  parameters  as  specified  by  (3.18). 
Pseudo-random  noise  is  added  to  the  data  in  the  manner  described  in  Example  6.2.  As 
seen  in  Figure  9,  the  Tikhonov  parametct  values  a  —  lO”'’  through  a  —  10“®  yield  points 
{\\K{ca)  —  d  II,  ||7?Cq||)  in  the  “corner”  of  the  L-curve.  For  Mx  =  16,  the  uniform  errors 
obtained  with  a  =  10~^,  a  =  10"*,  and  a  =  10"*°  are  reported  in  Table  4  with  correspond¬ 
ing  plots  of  the  true  and  approximate  parameters  shown  in  Figure  10.  As  indicated  by  the 
numerical  results,  the  “corner”  value  a  =  10“®  provides  a  good  choice  for  the  regularization 
parameter  whereas  a  =  10"*  causes  too  much  smoothing.  The  error  contributions  due  to 
the  smaller  singular  values  become  quite  apparent  at  cr  =  10"*°  thus  reiterating  the  Z-curve 
observation  that  this  Tikhonov  value  does  not  provide  enough  regularization  or  smoothing 
for  the  problem. 


a  =  10"*  a  =  10"*  a  =  10"*° 

!1pc/(0I1 

.7710  -  1  .4109  -  1  .6805  -  1 

Table  4.  Errors  on  the  Uniform  Grid  U  with  Mx  =  16  in  Example  6.4. 
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